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Abstract. We propose a mechanism to induce negative AC permittivity in the vicinity of a ferroelectric 
phase transition involved with spontaneous symmetry breaking. This mechanism makes use of responses 
at low frequency, yielding a high gain and a large phase delay, when the system jumps over the free- 
energy barrier with the aid of external fields. We illustrate the mechanism by analytically studying spin 
models with the Glauber-typed dynamics under periodic perturbations. Then, we show that the scenario 
is supported by numerical simulations of mean-field as well as two-dimensional spin systems. 

PACS. 42.65.Sf Dynamics of nonlinear optical systems; optical instabilities, optical chaos and complexity, 
and optical spatio-temporal dynamics - 78.20.Bh Theory, models, and numerical simulation (Optical 
properties of bulk materials) - 05.70.Jk Critical point phenomena 


1 Introduction 

When a many-body system is subjected to periodic driv¬ 
ing, e.g., via electromagnetic or acoustic waves, it has 
three different time scales: One is a microscopic time scale 
associated with thermal fluctuations, and we may regard 
this as a unit for measuring other time scales. Another 
time scale is related to internal collective dynamics, me¬ 
diated by interatomic coupling. The other is the period of 
the external driving. In the absence of periodic driving, 
the system will relax to equilibrium, and the fluctuation- 
dissipation relations provide a framework to relate the 
macroscopic relaxation to microscopic dynamics [I] . If ex¬ 
ternal driving is turned on, its period begins to com¬ 
pete with the other time scales. The interplay of these 
time scales has been investigated extensively in studies of 
stochastic resonance (see, e.g., Ref. [5] for a review). It is 
now well established that thermal fluctuations can play a 
constructive role to enhance sensitivity to the driving by 
modifying the internal time scale. 

When a spin system is perturbed by an oscillating field 
coupled to an order parameter, we have to observe the am¬ 
plitude and phase of the order parameter: In optics, for 
example, the optical response of an object to an incident 
beam is determined by the interference between the beam 
and secondary waves from the object, and their relative 
amplitudes and phase differences are important factors 
in the interference phenomenon [3]. By considering spin- 
spin correlation and its time scale in the generation of the 
secondary waves, we can obtain a more realistic picture 
of a real many-body system. Interestingly, the mean-field 
(MF) spin dynamics, the simplest description for collec¬ 


tive ordering, is similar to an overdamped oscillator, if 
written as an equation of motion for total magnetization 
as will be explained below. For this reason, the MF spin 
dynamics will serve as a starting point in our investiga¬ 
tion. Some progresses have been made through perturba¬ 
tive calculations whereby the existence of stochastic res¬ 
onance is predicted to the leading order [3115115] . One can 
also improve the prediction by taking higher-order terms 
into account [7|. However, we will point out a different 
kind of response that has been missing in these perturba¬ 
tive approaches. This motion is observed when the system 
hops from one free-energy minimum to another (see, e.g., 
instantons [8j for comparison). It shows a large phase de¬ 
lay relative to the incident wave, comparable to anoma¬ 
lous refraction. The frequency is nevertheless low. Such a 
large delay would be impossible at low frequency if the 
system was composed of weakly interacting simple har¬ 
monic oscillators. Moreover, we argue that if the system 
is effectively described as globally coupled, the response 
can be very strong compared with the input signal. It is 
known that the linear response theory can break down at 
low frequency wsm, and this failure has been related 
to an anomalously high gain puma. By considering 
another important factor shaping the response, i.e., the 
phase difference, we in this work will discuss anomalies in 
both the amplitude and phase, which we expect to open 
possibilities for interesting optical applications. 
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2 Model 

Some materials at a sufficiently low temperature exhibit a 
phenomenon called a polarization catastrophe, by which 
it acquires nonzero polarization P under a vanishingly 
small external electric field E E3- The consequence is 
that the static permittivity, defined as e = eo + P/E, di¬ 
verges at this point, where eo is the electric permittivity 
of free space. In other words, the system undergoes a fer¬ 
roelectric phase transition as temperature T varies. In the 
vicinity of the transition point T = T C , the behavior of the 
system can be phenomenologically described by the Lan¬ 
dau theory |151I16| . The central assumption is that the 
free-energy density can be expanded as a polynomial in 
the order parameter P: 

Fl = —EP + g Q + -cr(/3 c — P)P 2 + J9±P 4 + • • •, (1) 

where a and g 4 are positive constants, and f3 = (fc^T ) -1 
is the inverse temperature with the Boltzmann constant 
ks- In terms of /?, the critical point is now expressed as 
f3 c = {UbTc)- 1 . Minimizing F\, with respect to P requires 

0 = -E + u(/3 c - (3)P + g.iP 3 + .... ( 2 ) 

When E is absent, the theory predicts a continuous phase 
transition at f3 c of the MF universality class. Note that 
this is the universality class of three-dimensional quan¬ 
tum Ising ferromagnets and uniaxial dipolar Ising ferro- 
magnets [T7llTMT9| . 

Let us investigate the dynamical aspects by choosing 
the globally coupled Ising model as a concrete example of 
the Landau theory. We will discuss this model in magnetic 
terms as usual, but it should be understood as covering 
electric systems as well, if the magnetization order param¬ 
eter is substituted by polarization P. The Ising model has 
actually been used to describe ferroelectric materials such 
as NaNC >2 or (Glycine) 3 H 2 S 04 One might 

point out that polarization varies continuously in a crys¬ 
tal, making the naive Ising picture inappropriate. How¬ 
ever, the discreteness of an Ising spin plays only a minor 
role and it is the up-down symmetry that is more cru¬ 
cial [24j. Let us consider the energy function of the Ising 
model, 

= ( 3 ) 

i>j i 

where J > 0 represents the ferromagnetic coupling and Sj 
is a binary Ising spin variable that takes either of ±1 (see 
Ref. j25] for a connection to the antiferromagnetic case). 
The response of the spin system to time-dependent h is 
commonly studied by using the Glauber dynamics m, 
which shows qualitative agreements with experimental re¬ 
sults Under the Glauber dynamics, the transition 
probability from spin configuration g, to v , differing by 
a single spin, is defined as 

U '^ L 1 + exp[-/3([7 /J - U v )] ’ ^ 


where U v means the energy of spin configuration 77 = 
{si, S 2 , • • •, sn} consisting of N Ising spins, as defined 
in equation ([3]). By solving the master equation, one can 
readily derive the effective free-energy density functional 
as 

2 

jjl 

E = —-/JU 1 In cosh /?(to + ft) (5) 

with magnetization m = iV _1 JU ,s, (see Ref. [ 6 ] for the 
details). It is important that we have already taken the 
thermodynamic limit in equation ([5]) . so that this descrip¬ 
tion is free from any finite-size effects. Thermal noises 
are also taken into account in this description, because 
we work with the free-energy density functional. This dif¬ 
fers from a common approach through the time-dependent 
Ginzburg-Landau phenomenology in the form of a non¬ 
linear Langevin equation (see, e.g., Refs. mm, which 
contains a double-well potential and a noise term sepa¬ 
rately so that every observable should be averaged over 
the noise. We additionally note that only steady states 
are concerned throughout this work, so that initial tran¬ 
sients are assumed to have died out in our analytic and 
numerical considerations. In the context of ferroelectric- 
ity, the spin variable s* translates as an electric dipole 
moment, while m and h correspond to P and E, respec¬ 
tively. It is well known that the behavior of equation m is 
qualitatively the same as assumed in the Landau theory in 
equation ©■ The time evolution in the Glauber dynamics 
is formulated as 

(jrn f) T 

—— = — — — = — to + tanh/3(m + ft). ( 6 ) 

at dm 

This is purely relaxational dynamics with no inertia, as 
indicated by the absence of second-order derivatives, so 
it can be compared to an overdamped oscillator. In the 
limit of small h —> 0 , we may restrict ourselves to linear 
responses. Then, the relaxation time r is obtained as 


r = 


1 + /3(m* 2 



( 7 ) 


where to* > 0 denotes the magnitude of magnetization in 
equilibrium [B], which is nonzero only when (3 > j3 c = 1. 
If the free-energy functional for (3 > f3 c and ft, = 0 is 
approximated by 


E 


12 



*2 


) 


2 



( 8 ) 


we readily find that to* 2 ~ 3/3 _3 (/3 — /3 C ) via the Tay¬ 
lor expansion. This internal time scale r is an important 
quantity when a periodic external perturbation h(t) with 
frequency / = (27r) _1 o; is applied, because stochastic reso¬ 
nance takes place when the time-scale matching condition 
/ « t -1 is met to the leading order P1IT1[51IB1I25] . The 
minimum of E determines to* , and r is related to the cur¬ 
vature around it. We have one minimum when (3 < /3 C , 
while two minima appear when (3 > /3 C in the absence of 
ft. The important point is that each of them has a finite 
curvature, which allows the system to relax to an equilib¬ 
rium point within a finite time scale. At (3 = f3 c , on the 
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other hand, the local curvature vanishes around to* = 0, 
which is expressed by r —> oo in the linear-response theory 
(see Eq. 0). 

However, something is missing in this picture: When 
/3 > /?c? we can define an additional time scale t' in¬ 
volved with the swing between to* and —to*. This swing¬ 
ing motion becomes possible when the energy scale due 
to the periodic perturbation such as h(t) = ho cos uit over¬ 
comes the free-energy barrier in the middle. But the en¬ 
ergy scale needs not be large, because the barrier can be 
made as small as we want by adjusting the temperature 
close enough to the critical point. If the external perturba¬ 
tion barely overcomes the barrier, the swinging time scale 
t 1 can be much larger than the usual relaxation time r. 
It means the existence of a critical field strength h c above 
which the swing becomes possible, and t' diverges as h 
approaches h c from above. One can see the reason by ap¬ 
proximating equation © for small j3(m + h) as follows: 

~ (—to + tanh/3m) + /3h, (9) 

to the leading order in /3h < Cl. Let us assume that the 
external held h changes little while to transits from to* 
to —to*. The expression inside the parentheses is an odd 
function of m, so one may take the first sine Fourier com¬ 
ponent within an appropriate interval. Here, we simply ob¬ 
serve that it vanishes at ±to* and that its slope at to = 0 
equals (/3 — /3 c )m to approximate the functional shape as 


dm 

dt 


(/? - p c )m* 


sin 


r) +Ph. 
\m* / 


( 10 ) 


It is noteworthy that equation m can also be motivated 
by considering a gravity pendulum in the overdamped limit 
Enusmsa. The gravity pendulum should be inverted up¬ 
side down to describe one unstable fixed point (to = 0) 
and two stable fixed points (±to*). The swing time t' can 
be estimated from equation m as 


T 


/ 



dt 

dm 


dm 


2 m* 

/V^l’ 


( 11 ) 


where h c = TT~ 1 m*(l — f3 c //3) oc (T c — T) 3 / 2 . Clearly, the 
divergence of t' as h —>• hf has the same origin as the 
type-I intermittency [33] . 

Let us now take the temporal variations of the ex¬ 
ternal held into account. We will ignore the high-order 
terms of to in equation ©, because those terms are vis¬ 
ible only when |m| is sufficiently large, whereas we have 
argued above that the system spends most of the time 
in crossing over the barrier at small |m|. Based on this 
argument, we consider the following differential equation: 


dm 

dt 


(/3 — /3 c )m + (3ho cos cot, 


( 12 ) 


where 

* s MCt “ (yrs) ■ (14) 

Discarding other nonphysical solutions which blow up to 
infinity, we are actually relying on the existence of high- 
order terms, which seem to be neglected in equation (fT21) . 

Even this crude picture hints at some of interesting 
characteristics. First of all, the phase delay S = <p + tt 
lies between n and r so that the overall sign of m(t) 
is opposite to that of h(t). This simply means that h(t) 
should be positive to drive the system with m(t ) < 0 to 
the other side, and vice versa. However, such a large phase 
delay would not be seen if we regarded this system as an 
overdamped oscillator inside a well, whose phase delay 
is bounded from above by ^ |5]. Note also that we can 
practically quantify the phase delay 5 by calculating 

y" 

tan S = —. (15) 

X 

where we have used Fourier cosine and sine components 
defined as 


x' = ^ /cT TO (0 C0SU}t d M)> 
x" = ^ fo nsinu}t d M)- 


Let us recap this characteristic in terms of t'\ It takes 2 t' 
to make a full swing back and forth, and this should be 
shorter than the driving period / _1 to sustain the motion. 
Otherwise, the field will change the direction before the 
system overcomes the barrier. This can be understood as 
the time-scale matching condition for this mode. From this 
consideration, we derive the following inequality: 

2 t' < f~\ (17) 


At the same time, if the system follows the field too quickly, 
the phase difference will not be as large as suggested above. 
We are interested in a condition to have to < 0 for a suf¬ 
ficiently long time when h > 0, and vice versa. We al¬ 
ready know that the field will be positive for a half period 
(2/) _1 , and require that the duration of to < 0 should 
occupy more than a half of it, i.e., (4/) _1 . We thus arrive 
at the following inequality: 

7° dt 

(if)- 1 < / dm . (18) 

J-m * dm 


When h is close to h c , the right-hand side (RHS) is com¬ 
parable to t' (Eq. JH])), because the system spends most 
time in jumping over the free-energy barrier at m < 0. Af¬ 
ter some algebra, one can summarize these two conditions 
as follows: 

fc< f < 2 f c , (19) 


where /3 > /3 C . One obtains a physically acceptable solu¬ 
tion as 


m(t) 


P 

V (/?-/ 3 c ) 2 +^ 2 


ho cos (ut + <f> + 7r), 


(13) 


where f c = (87t) 1 \J (7t/3/i/to*) 2 — (j3 — /3 C ) 2 , and to* = 
a/3/3 _3 (/3 — /3 C ). In addition, noting that the calculation 
of t' in equation m approximates a sinusoidal driving to 
a constant held, we may identify h 2 in equation (fl9l) with 
a squared average, f Q n h g cos 2 (ut)d(ujt) = ho/2. 
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The discussion above tells us that this mode is possible 
at low frequency, if h — > /i+ (see Eqs. m and (Liill) ). In 
this respect, we can ignore w 2 in equation (fT51) compared 
with (/3 — (3 C ) 2 , which is kept constant, and the expres¬ 
sion inside the square brackets in equation (1131) can then 
be greater than 0(1). We argue the reason as follows: As 
soon as we bring the system from one minimum to the top 
of the free-energy barrier with h > h c , the system sponta¬ 
neously moves over to the other side. During this process, 
therefore, the traveling distance can be very large com¬ 
pared with h, because h c is an arbitrarily small quantity 
and the free-energy landscape is quite flat when /3 > ft c . 
It implies a possibility of high gain in the sense that the 
resulting change of the order parameter can be large with 
respect to the input field strength. Of course, it should be 
noted that this amplification factor is system-dependent, 
while the negative response to h(t) is universally expected. 

A low-dimensional system may be an example of such 
system dependence. In a finite-dimensional space, spatial 
variations of m over coordinates r also contribute to the 
free energy in such a way that 

•^Mr)] = J j + Vp [m(r)] j dr, (20) 



P 



where Vp can be identified with the RHS of equation ©. Fig. 1 . (a) Phase diagram of the MF kinetic Glauber-Ising 
By assuming relaxational dynamics to minimize equation CUlt nodel, driven by h = hocosut with ho = 10 -2 . We have nu- 
with respect to m( r), one derives a minimal description merically integrated equation © by using Heun’s method with 


with a Laplacian term added to equation © [34]: 

= V 2 m — to + tanh/?(ra + h). (21) 

This description assumes that the system lies deep enough 
inside the ordered phase to neglect critical-point fluctua¬ 
tions for the sake of self-consistency [23]. Put in the mo¬ 
mentum space, the equation transforms to 

dm(k,t) _ _(fc 2 -p i) m -p tanh/3(m + h), (22) 


time-step size At = 10 3 to obtain the solid lines. The phase 
delay 8 is estimated from equation m- The dotted lines are 
obtained by substituting h 2 = ho/2 into equation (1101) . (b) 
Time series of m responding oppositely to the external pertur¬ 
bation, when /3 = 1.04 in units of J -1 , and / = 4 x 10 -3 . We 
use the same ho and At as in panel (a). The solid line is ob¬ 
tained by integrating equation ©■ For comparison, the other 
curves show coscat and the approximate expression in equa¬ 
tion urn respectively. Note that the amplitude is much larger 
than ho = 10 -2 . One can also check from equation eni that 
x' is negative whereas x" is not, yielding 5 > 


where k = |fe|. We expect a characteristic length scale 
to exist in this system because /3 > /3 C . If that is the 
case, we may consider only this characteristic mode of 
nonzero k, which would correspond to the typical size 
of spin domains. Then, the amplification factor is esti¬ 
mated by the crude calculation above (see Eq. (fl3l) l as 
P/V (/3 — /3 C — k 2 ) 2 + w 2 , which cannot grow indefinitely 
due to the nonzero k, even if /3 is close to /3 C and ui is small. 
This argument suggests that the phenomenon reported 
here can be diminished in magnitude by domain dynamics 
that incorporates nontrivial spin fluctuations over space. 
For this reason, we speculate that the response of a low¬ 
dimensional system will not be as strong as in the globally 
coupled case. 


region inside the ordered phase at low frequency where S 
exceeds | (Fig. [TJt) . Note that equation (fl9l) qualitatively 
explains the shape of the region. Moreover, the response 
amplifies the input field with a high gain by an order of 
magnitude inside the region (Fig. [TJd). If further lowering 
/, we observe a discontinuous phase transition which is 
explained in the adiabatic approximation j7]. 

Let us check how the behavior is affected by low dimen¬ 
sionality. The two-dimensional (2D) Glauber-Ising model 
has an energy function 

U = -jY,s i s j -hY J S i , (23) 

(i,j) * 


3 Numerical results 

The characteristic large phase delay is indeed observed 
in our numerical integration of the MF dynamics, equa¬ 
tion ©• Plotting 6 in the (/?, /) plane, we find such a 


where ~}2(ij) runs over the nearest neighbors. It can also 
yield an inverted signal, but the amplitude is not so large 
as in the MF model (Fig. [2^,), as was explained by equa¬ 
tion (l22ll above. Whereas the region of S > -k/2 extends 
to high / in the MF Ising model (see Fig. [ljr), it is not 
detected at / > 2 x 10 4 in our numerical calculation of 
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Fig. 2. Response of 2D Glauber-Ising system driven by a peri¬ 
odic external field with amplitude ho and frequency / = 10 -3 , 
as a function of /3. These plots are obtained by simulating 
Glauber-typed Monte Carlo dynamics on an L x L square lat¬ 
tice with periodic boundary conditions. The vertical dotted 
lines represent the 2D critical point /3^ D ^ = ln(l + v / 2)/2 ~ 
0.44. (a) The negative dip in y', indicated by the arrow, is 
found above /3c 2D \ In this plot, the system size is chosen as 
L = 80. (b) We observe a region in which phase delay 8 lies 
outside [0, 7t/2] when ho = 0.08. The curves look like zigzags 
because we have sampled grid points on the (/3, /) plane with 
finite resolution. 



(b) L=80 . 
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the 2D model (Fig. [^Jo). Still, an important point is that 
the region of large phase delays survives when we increase 
the system size. 

As another example, we consider the 2D five-state clock 
model on an I x I square lattice with size N = L 2 . 

Its energy function is given as 

E = - J Si ■ Sj - h • Si, (24) 

{i,j) i 

where each spin s j = (cos 9j, sin 9j) at site j = 1 
has a discrete angle 9j = 2irnj/q with rij = 0 , 1 ,..., 4 . 
The magnetization of the system is defined as a vector 
sum m = TV -1 y . s j. A weak and slow driving in-plane 
field h is applied in a perpendicular direction to m with 
amplitude ho -C 1 and angular frequency u<l 1J7]. The 
model in equilibrium undergoes double phase transitions, 
only one of which at the lower T is accompanied by spon¬ 
taneous symmetry breaking. Although its critical proper¬ 
ties belong to a different universality class from that of 
the 2D Ising model :3S], Fig. [3] shows that the symmetry 
breaking can induce a response with a large phase delay, 



P 


Fig. 3. Response of 2D five-state clock model driven by a peri¬ 
odic external field with amplitude ho and frequency / = 10~ 3 
in the perpendicular direction to the magnetization. We plot 
x' as a function of /3 by simulating Glauber-typed dynamics on 
an 80 x 80 square lattice with periodic boundary conditions. 
In equilibrium, the model experiences spontaneous symmetry 
breaking around /3 « 1.1 [35j, above which the negative dip is 
located (the arrow). 


manifested by a negative dip in y/. as long as the field 
supplies enough driving force to overcome the small free- 
energy barrier. This is consistent with our understanding 
that the anomalous behavior is not specific to particu¬ 
lar model systems but related to general features of the 
symmetry-breaking phenomenon. 


4 Discussion and summary 

Although the existence of anomalous response with a large 
phase delay is precluded in the overdamped limit of the 
linear-response theory, we have found that it still remains 
possible as a non-perturbative mode. We have also ar¬ 
gued the characteristic features in the general terms of 
free-energy landscapes, meaning that it is quite a general 
mechanism involved with spontaneous symmetry breaking 
itself. The magnitude of the response depends on which 
specific system we are dealing with, and we expect the ef¬ 
fect to be better manifested if the system fits better into 
MF-like dynamics. In particular, if polarization P is oppo¬ 
site to the external field E with sufficiently large magni¬ 
tude, it implies that the AC permittivity e(w) = eo + P/E 
may become negative. Whether this effect is observable 
experimentally remains to be tested in a future study. 


S.K.B. was supported by Basic Science Research Program through 
the National Research Foundation of Korea (NRF) funded 
by the Ministry of Science, ICT and Future Planning (NRF- 
2014R1A1A1003304). B.J.K. was supported by the National 
Research Foundation of Korea (NRF) grant funded by the Ko¬ 
rea government (MSIP) (No. NRF-2014R1A2A2A01004919). 


References 

1. P. Hanggi, H. Thomas, Phys. Rep. 88, 207 (1982) 

2. L. Gammaitoni, P. Hanggi, P. Jung, F.Marchesoni, Rev. 
Mod. Phys. 70, 223 (1998) 



















6 


Seung Ki Baek et al.: Anomalous response in the vicinity of spontaneous symmetry breaking 


3. E. Hecht, Optics, 4th edn. (Addison Wesley, San Francisco, 

2002) 

4. K.T. Leung, Z. Neda, Phys. Lett. A 246, 505 (1998) 

5. B.J. Kim, P. Minnhagen, H.J. Kim, M.Y. Choi, G.S. Jeon, 
EPL 56, 333 (2001) 

6. S.K. Baek, B.J. Kim, Phys. Rev. E 86, 011132 (2012) 

7. B.J. Kim, H. Hong, J. Korean Phys. Soc. 52, 203 (2008) 

8. R. Rajaraman, Solitons and Instantons (North Holland, 
Amsterdam, 1987) 

9. J. Casado-Pascual, J. Gomez-Ordonez, M. Morillo, 

P. Hanggi, EPL 58, 342 (2002) 

10. L. Schmidt, R.R. Netz, EPL 98, 10014 (2012) 

11. S.K. Baek, F. Marchesoni, Phys. Rev. E 89, 022136 (2014) 

12. J. Casado-Pascual, J. Gomez-Ordonez, M. Morillo, 

P. Hanggi, Phys. Rev. Lett. 91, 210601 (2003) 

13. J. Casado-Pascual, C. Denk, J. Gomez-Ordonez, M. Mo¬ 
rillo, P. Hanggi, Phys. Rev. E 67, 036109 (2003) 

14. J. Casado-Pascual, J. Gomez-Ordonez, M. Morillo, 

P. Hanggi, Phys. Rev. E 68, 061104 (2003) 

15. C. Kittel, Introduction to Solid State Physics (John Wiley 
& Sons, New York, 1953) 

16. P. Chandra, P. Littlewood, in Physics of Ferroelectrics 
(Springer, Berlin, 2007), Vol. 105 of Topics in Applied 
Physics, pp. 69-116 

17. N.J. Als-Nielsen, L. Holmes, H. Guggenheim, Phys. Rev. 
Lett. 32, 610 (1974) 

18. N.J. Als-Nielsen, Phys. Rev. Lett. 37, 1161 (1976) 

19. G. Ahlers, A. Kornblit, H.J. Guggenheim, Phys. Rev. Lett. 
34, 1227 (1975) 

20. Y. Yamada, I. Shibuya, S. Hoshino, J. Phys. Soc. Jpn. 18, 
1594 (1963) 

21. Y. Yamada, Y. Fujii, I. Hatta, J. Phys. Soc. Jpn. 24, 1053 
(1968) 

22. R.E. Nettleton, Ferroelectrics 2, 77 (1971) 

23. T. Mitsui, E. Nakamura, M. Tokunaga, Ferroelectrics 5, 
185 (1973) 

24. J. Cardy, Scaling and Renormalization in Statistical 
Physics (Cambridge University Press, Cambridge, 2002) 

25. E. Vives, T. Castan, A. Planes, Am. J. Phys. 65, 907 
(1997) 

26. R.J. Glauber, J. Math. Phys. 4, 294 (1963) 

27. B.K. Chakrabarti, M. Acharyya, Rev. Mod. Phys. 71, 847 
(1999) 

28. S.G. Han, J. Um, B.J. Kim, Phys. Rev. E 86, 021119 

( 2012 ) 

29. B.A. Huberman, J.P. Crutchfield, N.H. Packard, Appl. 
Phys. Lett. 37, 750 (1980) 

30. E.G. Gwinn, R.M. Westervelt, Phys. Rev. Lett. 54, 1613 
(1985) 

31. C. Grebogi, E. Ott, J.A. Yorke, Phys. Rev. Lett. 56, 1011 
(1986) 

32. B.J. Kim, P. Minnhagen, P. Olsson, Phys. Rev. B 59, 
11506 (1999) 

33. Y. Pomeau, P. Manneville, Commun. Math. Phys. 74, 189 
(1980) 

34. P.L. Krapivskiy, S. Redner, E. Ben-Naim, A Kinetic View 
of Statistical Physics (Cambridge University Press, Cam¬ 
bridge, 2010) 

35. J.V. Jose, L.P. Kadanoff, S. Kirkpatrick, D.R. Nelson, 
Phys. Rev. B 16, 1217 (1977) 

36. S. Elitzur, R.B. Pearson. J. Shigemitsu, Phys. Rev. D 19, 
3698 (1979) 


37. H.J. Park, S.K. Baek, B.J. Kim, Phys. Rev. E 89, 032137 
(2014) 

38. S.K. Baek, H. Makela, P. Minnhagen, B.J. Kim, Phys. Rev. 
E 88, 012125 (2013) 

39. O. Borisenko, G. Cortese, R. Fiore, M. Gravina, A. Papa, 
Phys. Rev. E 83, 041120 (2011) 


